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Abstract 

We study the four-dimensional Zi random-plaquette lattice gauge theory as a model of topological 
quantum memory, the toric code in particular. In this model, the procedure of quantum error 
correction works properly in the ordered (Higgs) phase, and phase boundary between the ordered 
(Higgs) and disordered (confinement) phases gives the accuracy threshold of error correction. Using 
self-duality of the model in conjunction with the replica method, we show that this model has exactly 
the same mathematical structure as that of the two-dimensional random-bond Ising model, which 
has been studied very extensively. This observation enables us to derive a conjecture on the exact 
location of the multicritical point (accuracy threshold) of the model, p c — 0.889972 . . ., and leads to 
several nontrivial results including bounds on the accuracy threshold in three dimensions. 
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1 Introduction 



The lattice gauge theory has been studied mainly from the standpoint of particle physics. Many 
interesting results have been known for years for systems with spatially uniform couplings P^. The 
lattice gauge theory with randomness in couplings has not been a target of active studies mainly 
due to lack of physical motivation. However, it was recently proposed that the random lattice 
gauge model can be related closely with error-correcting processes of the toric code, an interesting 
example of quantum information storage (5] 0] . This observation serves as a strong motivation to 
investigate the random case systematically. 

In Refs. 0E1 E], the two-dimensional (2D) ±J random-bond Ising model (RBIM) and the 3D 
Z<i random-plaquette gauge model (RPGM) were studied in the context error corrections of the 
toric code. More precisely, the error-correcting processes of the bond-qubit model, where a qubit 
resides on each bond of the lattice, can be mapped to the RBIM. The plaquette-qubit model, where 
a qubit is located on each plaquette, is similarly related to the RPGM as far as error-corrections are 
concerned. 

In the toric code, we assume that the phase error and bit flip error may occur on each qubit due 
to decoherence and interactions with the environment. Fortunately, if the number of errors is not 
very large, we can determine the "error chains" or "error sheets" from the measured syndrome, or 
the positions where check operators give nontrivial signs (—1). We can then correct the errors by 
applying appropriate operators to the relevant error positions, and the encoded information, which 
is the simultaneous eigenstate of all check operators with trivial (unit) eigenvalues, is recovered. 
However, if the error rate becomes larger than a threshold (accuracy threshold), this procedure of 
error correction fails because nontrivial ambiguities arise in the error positions, and consequently 
the encoded information is lost. It is therefore very important to obtain the correct value of the 
accuracy threshold. 

The corresponding statistical models such as the RBIM or the RPGM have two parameters: the 
temperature T and the probability 1 — p for the interaction on a bond (plaquette) to be negative, 
which corresponds to the error rate in the context of toric code. The error-correcting properties of 
the toric code are related to the ordering behaviour of the RBIM/RPGM on a line in the p-T plane 
(phase diagram) . Error correction can be performed successfully in the magnetically ordered phase in 
the p-T plane of the RBIM (RPGM), which means stable storage of the encoded information, while in 
the disordered phase error correction fails and the encoded information is lost. The transition point 
(the multicritical point on the Nishimori line 0) is thus supposed to give the accuracy threshold 
of error correction. It is therefore crucial to know the correct location of the multicritical point of 
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these models for the design of reliable quantum storage, for example. 

These recent results have lead us to the study of the AD Z2 RPGM which can be related to 
the AD toric code. A part of the reasons to investigate the AD model, not 2D or 3D cases, is that 
the system with spatially uniform coupling is known to be self dual in AD 6 , which serves as a 
powerful tool of analyses. In addition, the model of plaquette toric code on the 3D cubic lattice 
under repeated measurement is also related to the AD RPGM. As has been mentioned above, the AD 
plaquette gauge model has a useful property of self duality, and we will analyze the AD RPGM using 
the duality technique. The formalism of duality transformation used in this paper is due to Wu and 
Wang 7 j . They proposed the method of duality transformation using Fourier transformation for the 
analysis of 2D non-random multi-component spin models. Recently it has been shown in Ref. jS] 
that the Wu-Wang duality is applicable to studying the 2D random spin model. We show that the 
duality analysis developed in Ref. [H] is also applicable to the study of the AD RPGM. Our main 
conclusion is that the mathematical structure determining the multicritical point of the AD RPGM 
between the ordered (Higgs) and the disordered (confinement) phases is closely related to that of the 
2D RBIM, whose property has been studied extensively for years. This suggests that the accuracy 
threshold of the AD plaquette toric code coincides with the 2D bond toric code. We also elucidate 
the correspondence between the AD plaquette qubit toric code and the AD RPGM. 

The 3D RPGM can be viewed as a model of the 3D toric code as well. We also study the phase 
structure of the 3D RPGM using the Wu-Wang duality and give useful bounds on the value of 
accuracy threshold of error correction. 

This paper is organized as follows. In Sec. 2 we review the toric codes for the 2D bond model and 
the AD plaquette model. We also elucidate the relation between the toric code and the random spin 
systems such as the RBIM and the RPGM. In Sec. 3 we explain the Wu-Wang duality technique 
using the 2D Ising model on the square lattice, which is one of the simplest self-dual models. In 
Sec. 4 self-duality of the AD plaquette gauge model without randomness is explained. In Sec. 5 
we incorporate randomness to the AD plaquette gauge model and analyze it using duality and the 
method of averaged Boltzmann factor in conjunction with the replica method proposed in Ref. [Hj- 
In Sec. 6 the 3D RPGM is investigated using its dual representation. In the dual space, the 3D 
RPGM can be transformed to an n-replicated spin system with spatially uniform couplings. Relation 
to the 3D toric code is also discussed. Section 7 is devoted to conclusion and discussions. 
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Figure 1: Qubits reside on bonds in the toric code (left). Possible error chains (right). 



2 Toric code in 2D and 4D 

For clarity of presentation, let us first review the toric code for quantum error correction and its 
relation to spin and gauge models. The toric code is a method to encode quantum information for 
its stable storage using topological nontriviality of the manifold on which qubits are located. 

It is useful to start with the 2D toric code to explain the basic idea |2J. In this case the toric 
code is denned on the 2D square lattice on the torus, and qubits reside on bonds (Fig. 1). We 
represent the two states of the qubit by 2-vector ( * J and ^ ° J . Any state of a specified bond 
can be expressed as a linear combination of these vectors. A quantum state of the whole system is 
given by the direct product of these states. The stored quantum information is chosen as a linear 
combination of such direct products which is a simultaneous eigenstate |W) of all check operators 
(to be defined below) with trivial eigenvalue 1. 

This system suffers from errors caused by decoherence and interactions with the environment, 
and the errors change the states of qubits on bonds. There are two types of errors; one is a bit-Hip 
error which flips the qubit to the other state, and the other is a phase error to change the phase of 
qubit. We express these errors by the operations of Pauli matrices to the original state. Existence of 
these errors on bond I is expressed as crf\^>) (bit-flip error) and erf | \& } (phase error). We can check 
if these errors have occurred by the operation of check operators. For this purpose, we define two 
kinds of check operators, 

Xi = (g)af, (1) 

iedl 

Zp = of, (2) 
ledP 

where Xi is denned on each site and Zp on each plaquette (shown on the left of Fig. 2). All check 
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check operators error correction 

Figure 2: Check operators and the inference procedure of error chain. Black circles denote syndrome. 

operators at any sites or plaquettes commute with each other, [X i} Zp] = 0, and hence simultaneous 
eigenstates of these operators exist. The eigenvalues of these check operators are ±1. We may there- 
fore choose the stored quantum state |\&) as a simultaneous eigenstate of all these check operators 
with the eigenvalue 1. 

Let us now consider the case that a chain of phase errors emerge. Bit-flip errors can be treated 
similarly on the dual lattice. We can detect the error chain by operating check operators to the 
stored quantum state. Phase error caused by of can be detected by the check operator Xi on the 
original lattice and a chain of bit-flip errors (due to af) by the operator Zp on the dual lattice. 
Since we prepare the original state whose eigenvalues of check operators are all positive (1), we can 
find the end positions of the error chain from the positions of a wrong sign (—1) of check operators, 
called syndrome. It should be remembered here that we can detect only both ends of the error 
chain, and we are asked to infer the actual configuration of the chain from its ends. In some cases 
our inference of error chain is successful, while in other cases we might infer a wrong error chain. 
An example of error correction procedure is shown in Fig. 2 (right). We can detect only the ends 
(black circles, syndrome) of a real error chain (E) as the sites where X^) — — |^) in the case of 
phase errors. An example of correction chain inferred from the syndrome is shown dotted (E 1 ). The 
example shown in Fig. 2 is a case of successful correction, where two chains E and E' are different 
but homologically equivalent. The reason will be explained later. 

In Fig. 3 are shown two examples of error correction. An example of successful correction is 
shown on the left. In this case, the real error chain (E, solid line) and correction chain (E 1 , dotted 
line) form a loop, which is homologically trivial (or can be shrunk to a point). Since the change of 
E to E' (or vice versa) is represented by the operators Z p on the plaquettes inside the loop E + E' 
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Figure 3: Successful and unsuccessful error corrections. 



and 1^) is invariant by their operation, the difference between E and E' is not essential, 

<gKI*> = (gKI*>- (3) 

leE ieE' 

An example on the right of Fig. 3 is an unsuccessful case. In this case, E and E' form a homologically 
nontrivial loop on the torus. It is impossible to change E to E' by operating Z p 's. Thus an error 
correction by using E' (dotted lines) yields a different (erroneous) result. Such an unsuccessful 
correction procedure often occurs when the error rate becomes large because many /long error chains 
may exist. It is expected from these examples that the accuracy threshold of error correction exists. 

Next we explain the relation between the 2D toric code and the 2D RBIM (Fig. 4) [3]. First let 
us consider the procedure of successful error correction when the real error chain E and correction 
chain E' constitute a homologically trivial loop. It is useful to consider the dual lattice. We assign a 
spin on each dual site and an interaction for each dual bond (Fig. 4, right). In doing so, we reverse 
spins inside the loop E + E' and assign antiferromagnetic (reversed sign) interactions to the bonds 
which correspond to the error chain E. By this procedure the 2D toric code can be seen as a (± J) 
RBIM. The error rate 1 — p of the toric code corresponds to the probability of antiferromagnetic 
interaction for each bond. The ratio of error probability to non-error probability (1— p) /p is identical 
with the edge Boltzmann factor of local interaction corresponding to the ratio of unfavourable and 
favourable spin alignment, e~ K je K as seen in Fig. 4. Thus the 2D RBIM under consideration lies 
on the Nishimori line defined by e~ 2K = (1 — p)/p 5 . 

The relation between the phase diagram of the 2D RBIM and accuracy threshold of the toric 
code is explained intuitively as follows 02]. In the magnetically ordered phase of the RBIM, error 
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Figure 4: Correspondence between the toric code and the RBIM. 

correction is successful since the islands of reversed spins as depicted on the right of Fig. 4 arc 
not extensively large in the ordered phase, which implies that the difference between E and E' is 
not significant. On the other hand, in the disordered phase, the error correction fails and encoded 
information is lost. Thus the critical probability of antiferromagnetic bond at the phase boundary 
along the Nishimori line (multicritical point) gives the accuracy threshold of error correction in the 
toric code. 

Now we generalize the toric code from the 2D bond-qubit system to the AD plaquette-qubit 
system. A motivation is that the accuracy threshold of the AD plaquette-qubit system can be 
estimated from the study of the corresponding AD gauge model which can be analyzed by duality, 
as shown in the following sections. For this purpose we prepare a AD hypercubic lattice and assign 
a qubit on each plaquette (Fig. 5). Bit- flip and phase errors occur on each plaquette. We define the 
check operator as follows, 

Xi = (g) o* P , (4) 

Zc = (g) o£, (5) 

Pedc 

where Xi is defined on each bond and Zc on each 3D cube on the AD lattice as shown in Fig. 5 
(right). We can check a phase error op by the operator Xi and a bit-flip error op by Zc- In the 
same way as in the 2D case, the phase error op is detected on the original lattice and the bit-flip 
error op on the dual lattice. Note that the plaquette model on the AD hypercubic lattice is self-dual 
in the sense that the dual system also carries its degrees of freedom on plaquettes as elucidated 
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Figure 5: AD toric code (I): Qubits errors and check operators are displayed. The 2D square lattice 
is a part of the AD hypercubic lattice. 

below. 

Let us consider the correction procedure of phase errors (Fig. 6). An error chain in the 2D 
bond toric code becomes an "error sheet" in the present AD case. We detect the error sheet using 
the check operator Xi. We can determine only the boundary of the error sheet from syndrome, 
which is similar to the 2D bond case. We must infer the error sheet itself from its boundary. In the 
correction procedure, the real error sheet E and the correction sheet E' form a connected surface. If 
the surface is closed or homologically trivial (middle in Fig. 6), error correction is successful because 
the operations of correcting operators on E and E' are equivalent to each other. On the other hand 
the correction procedure is not successful when the surface is homologically nontrivial (right in Fig. 
6). 

The correspondence of the AD plaquette toric code and the AD RPGM is similar to the 2D case 
3 . Here wc consider the case of successful error correction. Suppose that an error occurs with 
probability 1 — p on each plaquette. The probability Prob(£', E') that the error sheet E and the 
correction sheet E' are generated is given by (similarly to Fig. 4), 

Prob(£:,£') oc Y[exp(K P Up), (6) 
p 

where Up takes —1 on the surface C = E + E' and 1 elsewhere. Here Kp is defined by 

( = < j ( l ~P)/P> fOT P$ E ' (7) 

p/(l-p), for P EE. 
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Figure 6: The AD toric code (II): Error sheet and error correction procedure. The error sheet and 
correction sheet shown in the figure are hyperplanar surfaces. 

For each bond I on the AD lattice, Up must satisfy the constraint (because a connected surface C 
has no boundary), 

n ^ =i < ( g ) 

P-. leap 

or on the dual lattice, 

II U P' = !> ( 9 ) 

P': P'EdC 

for each dual cube C We introduce dual bond variables to solve this constraint, 

Up- (=U P )= \{ i v , (10) 
I'edP' 

where is an ordinary Ising variable which takes ±1. Using Kp'(= Kp) and Up>, we can rewrite 
the probability of Eq. © , 



Prob(£, E') cx Yl exp { K P , U P , (|,, ) } = exp | K ^ t p , U p , {fy ) j 



(11) 



where Tp> = Kp//K takes ±1. If we take the sum over This turns to the partition function of 
the AD RPGM, which we will study in Sec. 5. K is defined by e~ 2K = (1 — p)/p, the condition of 
the Nishimori line 0. 
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Figure 7: The 2D Ising model and its self-duality. Also illustrated is the derivation of modulo-2 
Kronecker delta of Eq. 1)21 [I for each plaquette. 



3 2D Ising model and Wu-Wang duality 

The next step is to develop a framework to investigate the mathematical structure of the models 
defined in Sec. 2. We consider the 2D Ising model on the 2D square lattice in order to illustrate 
the powerful technique of Wu-Wang duality. The Hamiltonian of the Ising model is 

H = -jJ2Si§j, (12) 

(ij) 

where Si is the normal Ising spin variable which takes ±1 and J is the uniform coupling constant. 
This Hamiltonian can be rewritten as, 

H = JJ2&, & = (13) 

l i£dl 

where I and i are bond and site, respectively, and Si is the modulo-2 spin variable which takes or 
1 in this case. The summation over i in the definition of is also defined by modulo 2. Note that 
a product of Si's has been changed to a sum of Si's. 
The partition function is, 

Z = E u(Z l )=eM-Ktii(Si)], (14) 

s ( =o,i i 

where K — (3J = J/kT. To take the dual of this model, we introduce a function obtained by 
regarding the modulo-2 bond variables as independent variables in Eq. 1)14(1 . 

~ z = E n<*)' «te) = cx P [-^6]. (15) 

£ ( =o,i I 



10 



Note that two functions Eqs. (|14|l and (|15J) do not coincide in their naive forms. To establish their 
equivalence, we must impose the following condition on the summation Eq. 1)15(1 . 

Y & = (mod 2) for all P, (16) 
leap 

where P is a plaquette on the 2D square lattice, because £j is composed of the modulo-2 sum of two 
Si's at the ends of /. This means that change of the value of spin Si at an arbitrary site on a 2D 
plaquette always causes modification two of the four £;'s in the plaquette, and l.h.s. of Eq. (|16[1 for 
any plaquette remains zero modulo 2. We denote the summation with this condition by a prime. 
For Eq. JEJ, 

Z{u) = const, x Y 

6=o,i i 

and we may then identify the two Z's in Eqs. I|14|) and 11711 . 

The dual representation of the partition function is derived as follows. The dual of a bond I (ID 
object) on the 2D square lattice is also a bond /' (Fig. 7). We define the dual bond variables r\y on 
the dual lattice, which also take the value or 1, and perform a discrete Fourier transformation of 
w(£z) defined by 



or conversely, 



*(vi f ) = -7= Y ex p( 7r *fwM6) = E ^("(E 

VZ 6=0,1 V 6=0,1 V Xiedi 

j= Y exp ™ PT 

v »),/=o,i \ \iedi 



5, m> (is) 



= ' - exp(^( y^Si] m > )u*M- Hi) I 



More explicitly, the edge Boltzmann factors are, 

u(0) = 1, 

u(l) = e~ K , 

u*(0) = ±(l + e - K ), 

u*(l) = -L(l- e - K ). (20) 

We use Eq. 119|l in Eq. I|14|) and take the sums over Si . Considering that a single site always appears 
at the ends of four bonds on the 2D square lattice, the exponential factors in Fourier transformation 
lead (modulo-2) to Kronecker deltas such as, 



Dd 2 



>€dP' J 



_ for all P', (21) 

\i' 

where P' is the dual plaquette on the 2D dual square lattice (Fig. 7). This condition is equivalent 
to Eq. (|16J) and we can replace the summation over "qy with the constrained one. 
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Figure 8: Bond variables and a plaquette. 



Now we can express the partition function using dual variables, 

Z(u*) = const, x Yi u *M- (22) 

77,/=0,l I' 

We have therefore obtained a direct relation between the two partition functions in Eqs. Ijl7(l and 

Z(u) = const, x Z(u*). (23) 

The transition point, if unique, is identified with the fixed point of duality transformation of the 
edge Boltzmann factor, 

u(0)=u*(0), u(l)=«*(l), (24) 
both of which lead to the same relation 1 

e~ Kc =V2-l. (25) 



4 4D Z2 lattice gauge model and Wu-Wang duality 

Let us move on to the AD Z2 gauge model. We can treat this model using duality in an analogous 
manner to the 2D Ising model. The Hamiltonian is 

H = jY,U p , U P =Y,& (~ II (26) 
p ledP ledP 

where P and I are the plaquette and bond on the AD hypercubic lattice, respectively (Fig. 8), and 
J is the uniform coupling constant. Here £j is the modulo-2 variable which takes or 1 and the 

1 If we start with the Hamiltonian Ea. 1121 . the transition point will be given by the relation e _2Kc = V2- 1. This 
difference is caused by the definitions of the coupling J in Eqs. 1121 and 1131 . which differ by factor 2. 
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Figure 9: Illustration of the condition Eq. Ij29|l : A 31? cube and Up's of six faces are shown. The 
operation of changing for one edge (shown in a thick line) always induces the changes of two Up's. 



sum over I is also defined by modulo 2. (£/ is the ordinary Ising spin variable which takes ±1.) The 
partition function is, 



Z= 11**^)' u(U P ) = e X p[-KU P (ti 

ii=o,i p 

To take the dual, we regard Up as the independent plaquette variables, 



(27) 



Z= n< U p)' u(U P ) = eM-KU P }. (28) 

u P =o,i P 

Functions (j27(l and (|28(l do not necessarily coincide unless we impose the following condition on the 
summation Eq. (|28|l . 

Y u p = ( mod 2 ) for a11 c i ( 29 ) 
Pedc 

where C is a 31? cube on the AD hypercubic lattice. This means that, if we change the value of 
at an arbitrary edge of the 3D cube, two of six Up's on the faces automatically change, and l.h.s. of 
Eq. (|29l) for any cube remains zero modulo 2 (Fig. 9). We denote the summation with this condition 
with a prime. For Eq. J5SJ, 

Z(u) = const, x Y Y[ u ( u p)- ( 30 ) 
u P =o,i P 

This expression is identified with Eq. I|27|l . 

Next we represent the partition function in the dual form. The dual of a plaquette P (2D object) 
on the AD hypercubic lattice is also a plaquette P'. We define the variables Vpi on dual plaquettes, 
which also take the value or 1, and perform a Fourier transformation of u(Up), 

u*(FpO = -L Yl exp{™V P Vp,)u(Up) = i Y expLi(j2b)vp') u ( U rl ( 31 ) 
V 2 Up=0yl V2 u p — o i \ \i e9P J ) 



or conversely, 



u(U P ) 



vz v P ,=o,i V \iedP ) / 



U*(Vp/). 



(32) 
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Figure 10: Illustration of the condition Eq. (|33|1 in the dual space: 3D cube in the dual space is 
localized in the x direction. 



We insert Eq. I|32[) into Eq. (|27|1 and take the sums over . Considering that a bond is always at the 
boundaries of six faces, the exponential factors in the Fourier transformation yield deltas (modulo-2) 
such as, 

" for all C", (33) 



Jmod 2 



\p'edc 



where C is the dual cube on the 4D dual hypercubic lattice (Fig. 10). This condition is equivalent 
to Eq. (|29[) and we can replace the summation over Vpi with the constrained one. 
We have thus obtained the partition function using the dual variables, 

Z{u*) = const, x Yi u *( v P')- (34) 

Vp/=0,1 P> 

The two partition functions in Eqs. Ij30|l and (|34|l are related as, 

Z(u) = const, x Z(u*). (35) 
From this relation, we can derive the transition point as in the 2D case, 

e- Ke =V2-l, (36) 
which is equivalent to the one of the 2D Ising model Eq. (|25|l . 



5 4D random plaquette gauge theory 

Next we explain the technique of dual transformation to deal with the RPGM following Ref. [Sj. 
They studied the 2D (±J) RBIM mainly using the Wu-Wang duality. Our method for analyzing 
the AD RPGM is essentially the same as in Ref. [H]- 
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From the motivation related to quantum error correction explained in Sec. 2, we incorporate 
randomness in the plaquette gauge model Eq. (|26|l . 



H = jJ2rpU P , 



(37) 



where rp is a plaquette-dependent quenched random variable. Here we assume that Tp takes 1 with 
probability p and —1 with 1 — p. 

We make use of the replica technique for averaging over randomness in rp. That is, we prepare 
n replica systems and take the n — > limit in the end to obtain physical quantities. After the 
configuration average over random plaquette variables, the partition function is determined by local 
non-random plaquette Boltzmann factors because the system then acquires spatial homogeneity. In 
the following we define the "averaged" plaquette Boltzmann factor for a given probability p. 

First we consider the simplest case of n = 1 . We define the plaquette Boltzmann matrix 2 as 







+ (1~P) 






X 


Xi 


p 


















Xi 


X 



(38) 



where k± 



xq means the averaged Boltzmann factor for the configuration Y[i Ui = 1 and x\ 



for Yli Ui = — 1. For a generic n, we can obtain the 2™ x 2™ matrix A n by a recursive procedure, 



P 



+ (l-p) 



K- 



K- 



= A n — 



(39) 



An-\ B n -1 

B n -i A n -i 

where B n is obtained from A n by replacing Xk —* Sfe+i- The factor Xk is defined by generalizing xq 
and X\ in Eq. (|38J) and corresponds to the averaged plaquette Boltzmann factor for the configuration 
Y[i Ui = 1 in n — k replicas and — 1 in k replicas. The explicit form of Xk is 

- k nJ + (1 -p) K+ k K- n - k . 



X k = pK + 

The averaged partition function depends on averaged plaquette Boltzmann factors Xk 

[Z n ] av = Z n (xo, x\, . . . , x n ), 



(40) 



(41) 



where [ ] av means random average. 

On the dual lattice, we can also define the dual averaged plaquette Boltzmann factor Xk* similarly. 
The Boltzmann factors for the original and dual systems are related by Eq. I|31(l . and accordingly 
the relations between Xk and x* k for n = 1 are given by replacing u and u* with x and x* respectively, 

V2xq = Xq + Xx, 



V2xl = 



Xq - X\. 



(42) 



2 In this section and Sec. 6 we define the Boltzmann factor by the normal Ising spin variable in Eq. 1261 instead 
of the modulo-2 variable 5; in order to follow the notation of Ref. 0. 
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For n = 2, 

2xq = (x + x\) + (.ti + x 2 ) — xq + 2x\ + x 2 , 
2x{ = ( x o - xt) + (xi - x 2 ) = x - x 2 , 

2x 2 = (x - xi) - (xi - x 2 ) = x a - 2xi + x 2 , (43) 

similarly to the 2D case 

In a similar way, we can obtain the explicit form of x* k for n replicas, 

2 n ' 2 x* 2m = ( K+ + K _)"- 2m ( K+ -^) 2m , 
2" /2 ^ m+ i = (2p-l)(« + + K -r- 2m - 1 ( K+ -«-) 2m+1 . (44) 
We can write the duality of the partition function using Xk and x* k , 

Z n (x ,xi, ...x n ) = Z h (xq,x*, . . .x*), (45) 

where we have neglected the trivial factor in Eq. (|35|) which is irrelevant to thermodynamic proper- 
ties. 

This procedure is completely the same as in the 2D RBIM. In the 2D RBIM, the multicritical 
point on the Nishimori line defined by 

e -2K = IzZ (46) 
P 

is supposed to give the lower bound of the probability p for ferromagnetic (ordered) phase. In Ref. 
[Hj, they conjectured that the phase boundary on the Nishimori line coincides with the crossing point 
of Eq. 146|l and a line defined by the relation for the averaged Boltzmann factor, 

x = x* Q . (47) 

They confirmed that this relation together with the condition of Eq. 146(1 leads to the exact multi- 
critical point at least in the case of small and infinite n. 

We write Eq. Q47JI explicitly with inverse temperature K and probability p, 



pe 



iK 



+ (l-p)e- nK = 2- n / 2 (e K + e- K ) n . (48) 



In conjunction with the Nishimori relation Eq. (|46|l . this yields, 

^ +1 + (l-p c )" +1 =2-"/ 2 , (49) 

where p c is the probability at the critical point on the Nishimori line (Fig. 11). If we expand Eq. 
(|49|l with respect to n and take the n — > limit, we obtain the relation of order 0(n 1 ) terms, 

-Pclogp c - (1 -p c )log(l -p c ) = (50) 
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p c =0.889972.. 



Figure 11: Phase diagram of the 2D ± J RBIM 



This equation gives p c — 0.889972 . . . numerically, which is in good agreement with the numerical 
studies in Ref. UJ. 

This argument for the 2D RBIM is also applicable to the AD RPGM. In this latter model the 
order parameter is defined by the Wilson loop operator, 

W[C] = \{lu (51) 

lec 

where C is an arbitrary closed loop. The random and thermal averaged value of W[C] obeys, 



ln{[W[C]] av , K ) cx { 



[length of C] (perimeter law) 

for ordered (Higgs) phase, 

(52) 

[area of C] (area law) 

for disordered (confinement) phase, 



where [ ] av ,K means random and thermal averaged quantity. Higgs (ordered) and confinement 
(disordered) phases correspond to ferromagnetic and paramagnetic phase, respectively. If we apply 
the previous discussion for the 2D RBIM to the AD RPGM, we automatically obtain the phase 
diagram for the AD RPGM from the one for the 2D RBIM by replacing the names of phases. The 
transition probability at the critical point on the Nishimori line is also given by Eq. (|50|l . We 
therefore conjecture that p c = 0.889972 . . . satisfying Eq. (|50|l is the exact accuracy threshold of the 
AD plaquette toric code corresponding to the AD RPGM. 
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6 3D RPGM and dual representation 



In this section we analyze the dual Boltzmann factor for the 3D RPGM and its representation by 
a spin model with ferromagnetic interactions. The discussion follows Ref. 8 where they mainly 
treated the 2D replicated ± J Ising model. 

We consider the ± J RPGM with the Hamiltonian 

H = -J^tpUp, (53) 
p 

on the 7i-replicated 3D cubic lattice. The probability distribution of t p is the same as in the previous 
section. The averaged Boltzmann factors x% and the dual averaged Boltzmann factors x* k are defined 
similarly to Sec. 5. It should be remembered, however, that the interpretation of dual Boltzmann 
factor is different in this case; In the 3D system the plaquette model is not self-dual (that is, the 
duality relation of the partition function in Eq. <|45l) does not hold), and the dual of a plaquette 
(2D object) is a bond (ID object) as illustrated in Fig. 12. If we consider the averaged Boltzmann 
factor Xk on a plaquette on the original lattice, its dual Boltzmann factor x k for a dual bond I* is 
defined for the configuration S*S* — 1 (S* is a dual spin variable and i,j € dl*) in n — k replicas 
and S*S* = —1 in the remaining k replicas. 

It is instructive to take the ratio of dual Boltzmann factors to Xq , 

(\ 2m— 1 
^+-^ =(2 J ,-l)tanh 2 '"- 1 ^, 
K + + ft_ J 

x* m /x* = tanh 2m if. (54) 

This ratio can be written in the simple form 

exp{K*(S* {1) + 5* (2) + . . . + S** (n) ) + K$S*WS*W . . . ST™} = Z n (K, K p ), (55) 

where tanhif = e~ 2K and 2p — 1 = e~ 2K p(= tanhK p ). S*^ is the Ising interaction factor in the 
fcth replica, 

5 *(fc) =s *(k) s *(k) (ij e dl*). (56) 

From the Boltzmann factor Ea. (|55[) . the dual model can be interpreted as a spin model which has 
ferromagnetic Ising interactions in each replica and an interaction between replicas. 

Next we take the K* — * oo (K p — * 0) limit which corresponds to the p — 1/2 case. In this case 
x 2m-i/ x o are zer0 f° r an y m an d S*^ S*^ S*^ . . . 5 , *(") is unity. Therefore a spin variable can be 
removed by the relation 5**(™) — S*^S*^ . . . 5 , *(™- 1 ) and the dual Boltzmann factor becomes 

Z n (K,K p = 0) 
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= Aexp{K*(S* w + 5* (2) + . . . + S* in - 1) + S* W S* (2) . . . S* (n - 1 '>)} 
= Z n ^(K,K p = K). 



(57) 



Thus the n-replicated system with K p — (p — 1/2) is equivalent to the (n — l)-replicated system 
with K — K p (on the Nishimori line). The Boltzmann factor Eq. I|55|) includes only ferromagnetic 
interactions and the Griffiths inequality |l()j holds, which leads to monotonicity of order parameters 
with respect to the change of K and K p (T and p). In particular the phase boundary between 
ferromagnetic and non-ferromagnetic phases is found to be a monotonic function of T and p (which 
are the parameters of the original model(=3£> RPGM)). This fact yields, 

T c "(p = 1/2) (= I^-^MCP)) < T c n (MCP), (58) 

where T™ is the critical temperature for the n-replicated system and MCP means the multicritical 
point. We rewrite this inequality using probability p, 

p^OMCP) >p£(MCP), (59) 

where p™ is the probability of the multicritical point for the n-replicated system. Let us suppose 
that this inequality holds in the n — > limit. If we consider the n = 1 system on the Nishimori line 
[K* = K*), the dual Boltzmann factor becomes, 

Z X {K,K V = K) = Aexp{2irS"* (1) }. (60) 

This is the simple ferromagnetic Ising Boltzmann factor and we can make use of the knowledge 
of the 3D ferromagnetic Ising model for the analysis. The critical point of the 31? non-random 
ferromagnetic Ising model is estimated as 2K* ~ 0.22165 . . . numerically Using the condition 
(K = K*) and the dual relation (2p - 1 = e~ 2K p), we obtain the value pJ(MCP) ~ 0.9006. Finally 
we obtain the following inequality, 

p°(MCP) > Pc(MCP) ~ 0.9006. (61) 

This inequality gives a lower bound of probability at the multicritical point for the 31? RPGM under 
the assumption that the inequality holds in the limit n — > 0. 

The 3D RPGM is related to the toric code as well. We consider the toric code on the 3D cubic 
lattice and locate a qubit on each bond. If we choose check operators for this system as shown in 
Fig. 12 (bond case), the accuracy threshold of phase error correction is given by the probability at 
phase boundary (on the Nishimori line) in the 3D RPGM. The accuracy threshold of bit-flip error 
correction is given by the phase boundary in the 3D RBIM. On the other hand, if we locate a qubit 



19 




"71 



X s Z P 

(bond qubit model) 



(plaquette qubit model) 



Figure 12: Duality of the 3D cubic lattice and check operators: X s (Xi) is given by the product of 
Pauli matrix a x for each bond (plaquette). Zp(Zc) is given by the product of a z in the same way. 



on each plaquette and choose the check operators like in Fig. 12 (plaquette case), the accuracy 
threshold of bit-flip error correction is given by the probability at phase boundary (on the Nishimori 
line) in the 3D RPGM. Similarly, the accuracy threshold of phase error correction is given by the 
critical point in the 3D RBIM. From this correspondence, we easily find out that the procedure of 
phase error correction in the 3D bond qubit system and that of bit-flip error correction in the 3D 
plaquette qubit system are related by duality. In the same way, the bit-flip error correction in the 
3D bond system and phase error correction in the 3D plaquette system are dual to each other. 



7 Conclusion and discussion 

We have discussed the phase structure of the 4D RPGM using the method of duality and averaged 
Boltzmann factor. We have derived an equation to determine the location of the multicritical point 
(the accuracy threshold in the context of quantum error correction) , which is expected to be exact 
(but is a conjecture, rigorously speaking). The structure of duality relations and the resulting 
equation for the multicritical point coincide precisely with the corresponding ones for the 2D RBIM. 
We could therefore generalize the generic duality relation originally proposed by Wegner |B] for non- 
random systems to the random case. In particular, the present argument can be extended to the 
model with n-dimensional objects in the 2n-dimensional lattice. The same result will be able to be 



20 



Nishimori Line 



Disordered 

p=0.889972 



-\ 1- 




0.5 0.6 0.7 0.8 0.9 1.0 



4D RPGM 



v Nishimori Line 



(B) 



p=0.9006... (D) j 



Disordered 



^Ordered 



2.0 



M7313...(F) 

1.0 



p=0.889972... (C)T 



0.5 0.6 0.7 0.8 0.9 ♦ 1.0 



V 



o.o 



0.9707... (E) 



3D RPGM 



Figure 13: Phase diagram of the AD and the 31? ± J RPGM. Point (A): Critical point of non-random 
system obtained exactly from the duality of non-random model. Line (B): Phase boundary of the AD 
RPGM. Line (C): Lower bound of p from the analysis of the AD RPGM. Line (D): Lower bound of 
p from the analysis of the dual model of the 3D RPGM in Sec .6. Point (E): Critical point at T = 
by the numerical study in Ref. 0|. Point (F): Critical point of the non-random model obtained by 
the numerical study of the 3D Ising model and by the duality between the 3D Ising model and the 
3D plaquette gauge model. 



derived for these models. 

We have elucidated the equivalence between the AD RPGM and the AD plaquette qubit toric 
code (and similarly for the 3D plaquette system). The phase boundary between the ordered and 
disordered phases on the Nishimori line gives the accuracy threshold of the toric code. From our 
study of the AD RPGM, the accuracy threshold of the toric code with qubits on AD plaquettes is 
expected to be given as p = 0.889972 . . ., the same result as that of the 2D bond qubit system. 

We may also expect that the phase boundary of the AD RPGM gives a limit to the boundary of 
the ordered phase for the 3D RPGM (Fig. 13) because the reduction of dimension always enhances 
disorder. Of course the phase boundary of the AD RPGM gives the lower bound of probability of 
ferromagnetic phase for the 3D RPGM (p = 0.889972 . . .). The lower bound was also estimated from 
the 3D Ising model using duality. The Griffiths inequality gives the lower bound p — 0.9006 ... on 
the assumption that the inequality holds for n — > limit. This is a better bound than the one from 
the study of the AD RPGM (p = 0.889972 . . .). The transition point at T = (which is expected 
as the probability at the multicritical point) is estimated at p = 0.9707 ± 0.0002 in Ref. : 4, by a 
numerical study, which is consistent with our inequality. 
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As for the AD plaquette toric code, a local error correction procedure is also proposed in Ref. 
[3]. The method is to be contrasted with that in the present paper where we use global information 
to infer the error positions because one needs the information of the boundary of error sheet. The 
difference between these correction procedures is as follows. Global procedure needs processes of fast 
classical computation to infer the real error sheet from its boundary, while the local procedure does 
not require such a classical computation because the local procedure needs only the local information 
of defects to remove all error sheets. On the other hand, the global procedure involves only a single 
observation to infer the error sheets and to correct errors. (Of course the observational errors occur 
in actual procedure and we must take them into consideration.) The local procedure requires many- 
step observations to erase all error sheets completely: The single correction step reduces the area 
of error sheet, but not to wipe them out at once. At the same time, we must notice that new 
error sheets may be generated during this process. If the rate of error generation is larger than the 
error reduction rate of the correction procedure, the error correction will fail because homologically 
nontrivial error sheet will be formed by the accumulation of error sheets, which affects the encoded 
information. (Remember that the homologically nontrivial sheet is formed by real error sheets and 
correction sheets in the global procedure.) Therefore the accuracy threshold also exists in the local 
procedure. In Ref. 3 an upper bound of accuracy threshold of the error correction by the local 
procedure is estimated from the argument mentioned above, p c ~ 1 — (4.8 x 1CP 4 ). As pointed out 
in Ref. [2], the local correction procedure cannot be applied to the bond-qubit system. We need at 
least four dimensions in order to perform the local correction procedure both for phase and bit-flip 
errors. 
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